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A method is developed which speeds up averaging in quantum simulations where minus signs 
cause difficulties. A Langevin equation method in conjunction with a replication algorithm is used 
enabling one to average over a continuously varying complex number. Instead of ensemble averaging 
this number directly, the phase of the complex number is followed over time. The method is 
illustrated in some simple cases where the answers obtained can be compared to exact results, 
and also compared to conventional averaging procedures which converge orders of magnitude slower 
than this method. Limitations of this method are also described. 



I. INTRODUCTION 

Simulation of quantum many body systems have been greatly hindered by the "sign problem" . Quantum systems 
with many degrees of freedom are often simulated using stochastic methods. For a variety of problems, usually 
involving fermions, an average over many negative and positive numbers must be taken. Sometimes the sign does 
not effect results but in many interesting cases the resulting average is very small because of cancellations, so 
that in order to obtain good statistics, the number of realizations that must be taken is exponential in the inverse 
temperature, making it infeasible to get information about ground state properties of a system. 

The main idea of this work can be understood by considering the following example. It will be seen below that it 
describes the quantum mechanics of n bosons restricted to one site. Consider the stochastic equation for z(t) which 
is a complex function of time: 

z = if(t)z (1) 

Here f(t) is real gaussian noise with correlation function {f(t)f(t')) — vS(t — t'), where the angled brackets denote 
an average over the noise f(t). Choosing z(0) = 1 the above equation is easily solved giving 

z = e m (2) 

where 

9(t) = f f{r)dr (3) 
Jo 

Here 6(t) describes a random walk in time, and therefore the the probability distribution of 9 at time t, P(6,t), is a 
gaussian distribution with zero mean and variance a 1 = vt. Below we will see that the ground state energy of n 
bosons is simply obtained from the asymptotic behavior of (z(t) n ), which is solved for as follows: 

(z(t) n ) = (e lne ) = f P(6,t)e me d6 = e -« 2 "*/2_ (4) 



Now let us turn to the problem of how, without the aid of this analytical solution, one could obtain (z(t) n ) 
numerically. Because \z(t)\ = 1, to obtain the correct average would require over e n vt independent runs of eqn. (|l|) 
to obtain adequate statistics. This can be very large. 
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The alternative to this direct averaging procedure is to obtain the probability distribution P(9, t). In this case it is 
much better to keep track of the total angle 9, not 9 mod (2tt). Keeping track of the total angle is done as eqn. (|l|) 
is evolved. The increment in 9 found in a single time step is added to the previous value of 9 and a histogram of 9 at 
a given time can be obtained by running eqn. (]l|) over many realizations of /. 

Does analytically continuing 9 and then finding P(9,t) help in obtaining (z{t) n )l In this case and in the less trivial 
cases considered below, it appears to reduce the computation time from something exponential to algebraic. Of course 
assumptions must be made in order to obtain such a reduction, the most important one being that P{9, t) is a smooth 
function of 9. There are also problems encountered in deciding what function to use in fitting P, as will be discussed 
later. In the case of interest here, fitting In P(9,t) to a quadratic function of 9 gives accurate results. In fact the 
ground state energies should be obtained to within an error A -1 / 2 where N is the number of independent runs. 

Having motivated the method that will be considered, the Langevin equation that simulates quantum many body 
systems will be described and numerical results will be given. 



II. SIMULATION METHOD 



A system of fermions or bosons can be simulated by considering a stochastic equation which is a generalization of 
eqn. (|l|). A field (j)(r,t) is governed by the equation 
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(r,i) = /(r,i)0+— V 2 < 



(5) 



where (/(r, t)/(r', £)) = S(t — t')v{v — r'). This describes a system of bosons of mass m interacting with potential 

n 

U = -]T w(ri - 1^-^(0) (6) 
evolving in imaginary time. The wave function at time t, is 

*(n,...,r„) = (0(n)...0(r n )>. (7) 



To simulate n fermions, n fields 0i , </>2 , . . . , </> ra are evolved starting from different initial conditions but with the 
same /(r, t). The wave function at time t in this case can be expressed in terms of the average of a Slater determinant: 



*(ri,...,r„) = 



0i (ri) ■■■ 0i (r„) 



n (ri) ... n (r n ) 



(8) 



Note that the potential in eqn. (^|) is the negative of the correlation function. Therefore to simulate repulsive 
potentials requires an f(r,t) that is complex [gj. This method is closely related to a simulation method for the density 
matrix ||,|| and also to projector Monte-Carlo ||. 

In order to obtain the ground state energy for quantum systems such as the ones described above requires deter- 
mining the asymptotic exponential decay of the wave function. In principal this can be done by averaging over many 
realizations of /, however there are two reasons why this is impractical. 

The first is a problem of importance sampling. As an example, consider determining the ground state energy E(n) 
of n bosons. This is obtained by calculating < <j> n (r, t) > and fitting this to an exponential exp(E(n)t) for long times. 
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The problem is that E(l) ^ E(n)/n. As a result the typical behavior of <p n (r,t) and the average behavior are very 
different. What dominates the average are some very rare realizations of /. Another way of describing such behavior 
is in terms of the notion of intermittency [||. In order to get the correct average, one must use importance sampling. 
There are two ways to do this. One is to use Monte-Carlo keeping track of all time paths ||. The second is to use a 
replication algorithm (t]]|]. In order to keep track of angles, it is much more convenient to use the second method. 

A large number of copies of systems described by eqn. (|^) are run in parallel. The number of configurations that 
are replicated is proportional to their weight. For example, if we wish to measure (|</>| n+1 )/(|0|™) at some site r, we 
should choose a weight proportional to \cp\ n . The proportionality factor is chosen so that the number of copies of the 
system stays almost constant. This method is similar to others described in detail in refs. [Q and ||. 

The second problem is one encountered when simulating repulsive bosons and fermions. Here non-positivity leads 
to the problem discussed in the introduction, that an exponentially large number of runs must be taken to obtain 
adequate statistics. The phase method proposed here is designed to reduce the number of runs. This will be illustrated 
numerically for a system of repulsive bosons. One should note that although it is possible to choose a method for this 
problem where there are no minus signs this is an instructive example. It will be seen that the method used 

generalizes to fermions. 



III. OBTAINING THE PHASE 

Consider n bosons on a one dimensional lattice with L sites interacting via an on-site repulsive interaction U. This 
has been the focus of recent investigation using "world line" Monte Carlo [p|JlO|] , which is clearly better suited 
to this problem than the method here. We use the method above only to illustrate techniques to cope with negative 
signs. In obtaining ground state energies, any component of |\? > can be measured. As will be seen below, in order 
to make the phase method work it is important to make a judicious choice for this. We start by choosing to measure 
\E'(xi,X2, ■ ■ ■ ,x n ) where the x^'s are all different sites. Later we will see why using the same sites does not work as 
well. The prescription for computing ^ is given by eqn. (Q). 

The replication method outlined above was the case where all the Xj's are the same. To generalize the discussion 
above, one just chooses a weight w = \4>{xi)(f){x2) • • • 4>(x n )\. The ground state energy of n + 1 particles can be easily 
computed by knowing 

< 4>(xi) . . . cj>(x n+ i) > 



= <e l <V(x n+1 )| > w = [p(Q > t)e i@ dQ (9) 



< \4>{xi) . . .<f>(x n )\ > 

where the last average is done with respect to the weight w, and 9 is the total phase angle of 4>{x\) . . .<f>(x n -\-i)- 
P(Q, t) ee (\<j) n+1 \6(Q - 6{t))) w is closely related to P{6, t) described below eqn. ||. 

To test this out numerically, the number of lattice sites was chosen to be 8, and we started out by using a small 
value for U = 0.5. An average of 4497 copies were replicated. From this data P(Q,t) was obtained with good 
statistics. In order to obtain an improved estimate of the tails of this distribution, the replication algorithm was 
altered to weight large O, that is by letting w — > w exp(consi.O). By changing the constant in the exponent, P(0,t) 
was probed for different regions of O. By combining these, one obtains a good estimate for P(Q,t) over nine orders 
of magnitude. An example of lnP(©,£) is shown in fig. [3] for n — 8. The solid line is a quadratic fit to the data. 
By taking the fourier transform of this for different times, as prescribed by eqn. (||), the ground state energy can be 
obtained. This is shown in fig. ^ by the solid triangles. The open squares show the exact results for comparison. 
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FIG. 1. The histogram P(8,t) at t = 5. The histogram was obtained by using the replication algorithm in conjunction with 
eqn. (||) . The algorithm was run by giving favorable weights to large angles enabling good statistics for large values of the 
angle to be obtained. The solid line is a quadratic fit to the data. 
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FIG. 2. The ground state energy E, as a function of number of particles n, for repulsive bosons with interaction energy of 
0.5. The solid triangles are the results obtained using P(Q,t). The open squares with the dashed line going through them are 
the exact results. 



Without using P(Q,t), one can compute the weighted average of \<j){x n+ i)\ exp(z6) directly. For comparison the 
average of this, over the same copies, is shown fig. ^|. It is clear that results obtained by direct averaging are far 
inferior to those obtained using P(Q,t). 
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FIG. 3. The result of averaging the field over the same runs as was used to obtain fig. [l| An exponential decay is not 
discernable and many more averages would be needed to see it. The real part is plotted here. 



IV. IMPROVEMENTS TO THE METHOD 

The method as it stands has difficulties for large U. The histogram of 8 has too large a variance. This can be 
understood by considering a variant of eqn. (|l|) in which there are two independent random processes with different 
/'s, (/i(t)/i(i')) = viS(t — t'), and (f2{t)f 2(f)) — v^Sit—t'). The equations corresponding to each are Zj = ifi(t)zi. 
and therefore Zi — e 10 ^*', % = 1, 2 . Suppose we are interested in numerically computing the sum 

{Aiz? + A 2 z%) = A ie - n2vit / 2 + A 2 e- n2v2t/2 (10) 

This average will be dominated for long times by the Zj with the smaller of the ViS, say v\. Under what circumstances 
will the numerical method of finding the histogram for the phase angle of the above sum give a sensible answer? Eqn. 
( |l0| ) can be thought of as the sum of two randomly rotating vectors with magnitudes \Ai\. From this observation, it 
is easy to see that for long times, the variance of the phase angle will be the same as that for z™ if \A\\ > \A2\. For 
l^-i I < 1^4-2 1 a, broader distribution will be obtained. Unless very careful fitting to a non-quadratic expression is done, 
an erroneously high energy will be obtained. 

This example shows that very careful fitting is necessary if the wave function has too high an amplitude of excited 
states. One way to reduce this effect is to choose to measure \& on different sites as was done in the numerical example 
above. Another way to improve this method is to filter (j> in time. It is easily shown that applying a linear filter to 
(j>(t), j(t) = (f)(t)*g(t) does not change the ground state energy if the filtering function g(t) decays sufficiently rapidly. 
For example one can evolve the equation 

7(r,i) = -c 7 (r,t) + 0(r,i) (11) 

simultaneously with evolving eqn. (^) and measure 7™ instead of 0". This method works quite well. The energies 
for an eight site lattice for the case U — 4 are shown in fig. ^. The energies are in good agreement with the exact 
results up to n — 6 after which it becomes apparent that the higher energy states start to dominate the answer. 
The method will presumably improve if in addition to filtering, different sites are used as in the example above. 
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FIG. 4. The ground state energy E, as a function of number of particles n, for repulsive bosons with interaction energy of 
4.0, obtained using linear filtering of the field (solid triangles), and the exact energy (open squares with the dashed line going 
through them). 




FIG. 5. A log polar plot of the Slater determinant denned in eqn. (g) as a function of time, for three particles on a lattice 
of 4 sites with a repulsive potential of 2.5. 

Lastly, the case of fermions has been considered. The determinant formulation of this problem can be used to 
obtain a continuous phase as a function of time. Fig. (|^) shows the complex determinant dona log polar plot, that 
is with radius — ln|d(i)| and angle arg(d(t)), for 3 particles on a 4 site chain, with U = 2.5. The phase randomly 
rotates, similar to the case of bosons. 

In conclusion, a new method for dealing with cancellations in quantum systems has been developed and tried out 
in some simple cases. Further theoretical development of the determination of P(0, t) could greatly improve fitting. 
Recent theoretical work jl2| relating the phase used here to the Berry phase for smooth paths should be further 
explored in the context of this simulation method. It would also be interesting to consider more difficult problems 
such as the two dimensional Hubbard model with repulsive interactions using this method. 
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